ANL-HEP-PR-01-025 
|hep-th/0612158| 



Numerical Program for Computing <I> 3 

Amplitudes 

Gordon Chalmers 

UCLA Physics Department 
405 Hilgard Ave 
Knudsen Hall 
LA, CA 90095-1547 
e-mail: chalmers@physics.ucla.edu 



Abstract 



A computing program in Matlab is given that computes amplitudes in 
scalar (j) 3 theory. The program is partitioned into several parts and a simple 
guide is given for its use. 



The scattering amplitudes in phi 3 are computationally difficult at high 
loop orders due to the number of diagrams involved. Their number grows as 
(n — 1)! at n-point and at tree level. An automated program is given here 
that computes an arbitrary tree amplitude and then sews them together to 
obtain the quantum amplitudes. The program is written in Matlab and is 
given in several pieces. 

The input parameters for the n-point amplitudes are given in the subpro- 
gram phi3compute.txt. 

The subprogram NMatrix.m must be loaded with the number of external 
lines ranging from npointlower to npoint. These numbers span the numbers 
of external lines which are used in the trees that are sewn together into the 
loop amplitudes. 

The subprograms NodalComplex.m, treecall.m, innerproduct.m, and de- 
nomcall.m are saved as files in the directory to which the program can call. 

The program RainbowCompute.m is used to call the MonteCarlo simu- 
lation that computes the multi loop amplitude. 

The input variables are self-explanatory in the subpart Phi3Compute.txt. 
They are: kmomol, kmomor delimit the four-point kinematics to the four- 
point amplitude, gnumber is the power of the coupling constant, d is the 
dimension, kstep and sample are the step size to the partitions of momenta in 
the loop and sample is the number of samples of the internal momenta taken. 
In the program RainbowCompute.m the variables are MonteRunSet denoting 
the number of samples, kmax is the discretized momentum, sampleX is the 
number of external momenta data points. The output is delivered in the 
array RainbowCompute and there is a subroutine to plot the output. 

I typically call all of the programs and then paste in the subprogram 
RainbowCompute.m. The time to obtain the output really varies depending 
on the parameters and the number of internal loops. 

Background algorithms to the computations are obtained in [TJ, [2]. The 
program is contained in the tex source following the bibliography. 
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mmmmmmmmmmmmmmmmmsmsmmmmmmo 
mmmmmmmmmmmwmmwmmmmmmmmmmi 
%% 

%% Phi3Compute . txt 

%% 

% momentum conservation at every node 

% check the redundance of the phi3 graphs 

%% this is the header to phi3 computation 

%% be sure to preload the N matrix first with NMatrix ( npoint ) program 

gnumber = 6; 
ppoint =2; 
qpoint =2; 
d = 4; 

clear kmomol ; 
clear kmomor ; 
kstep =2; 
sample = 100; 

for i = 1:4 

kloffsct ( i ) = .02323232323; 

kroffset(i) = .01515151515; 
end 

kmomol (2 ,4 ,100) = — .004* kstep * sample / 2 ; 
kmomor (2 ,4 ,100) = — .004* kstep * sample / 2; 

for samplcx = l:samplc 
for i=l:d 

for j=l:ppoint 
kmomol ( j , i , samp lex )=kstep *rand+klof f se t ( j ) ; 

end 
end 
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for i=l:d 

kmomol (1 , i , samplex) = kloffset ( i ) + .004* kstep *samplex ; 

kmomol (2 , i , samplex) =(1— kmomol (1 , i , samplex) *kmomol ( 1 , i , samplex ) ) " ( 1/2) ; 
kmomor (1 , i , samplex ) = kroffset(i)+.004*kstep * samplex ; 

kmomor (2 , i , samplex ) =(1— kmomor ( 1 , i , samplex ) * kmomor ( 1 , i , samplex ) ) " (1/2) — 
kmomol ( 1 , i , samplex )— kmomor ( 1 , i , samplex )— kmomor ( 2 , i , samplex ) ; 

end 
end 

clear RainbowComputeOne ; 

RainbowComputeOnc=RainbowCompute( sample , gnumber , ppoint , qpoint , kmomol , kmomor ) 

for j =1:100 

tcstx ( j )=j ; 

end 
figure 

plot (tcstx , RainbowComputeOne ) 

nmmmmmmmmmmmmmwmmwmmmmmmmmmmmm 
%% 

%% RainbowCompute.m 

%% 
%% 

function [ RainbowCompute] = RainbowCompute( sample , gnumber , ppoint , qpoint , kmomol , 
kmomor ) 

% had to preload the Nmatrix containing the variable N 
% docsnt contain the asterisked data immediately below 
% function docsnt work so paste in the information 
% gnumber = 10; 
% ppoint=2; 
% qpoint=2; 

[ lincconfig , NodcNumbcrPcrm] = NodalComplex (gnumber , ppoint , qpoint ) ; 
load Nmatrix 
sampleX = 100; 
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%% kmax is the max value of the momenta individual components in the 
integral 

kmax = .01; 

% need an offset for the momentum to isolate uv from ir regime 

samplex=samplc ; 
m=l; 

d = 4; 

% Two Nodes: I nt er nalL ines+ppoint — 2+Int er nalL ines+qpoint — 2=gnumber 
InternalLinesMax = (gnumber+4— ppoint— qpoint ) /2; 
InternalLinesMin = 2; 

NodeNumbcrMax = (gnumber— ppoint —qpoint +2) / 2 ; 
%% the number of random samples in the integral 
RainbowCompute( sample ) =0; 
RainbowComputetwo=0; 

%% number of external momenta data points 
samplcX = 100; 

%% number of integral momenta data points 
MontcRunSet=100; 
kmeasure=kmax/MontcRunSet ; 

WMWMWMSWWo symmetrize in the momenta %%%%%%%9®Wo 

%% 

plr ( ppoint+qpoint ,d,100)=0; 
testvec (ppoint ,d,100)=0; 

%% repeat many times to montc carlo simulate the integrals 
for samplex = l:sampleX 
%% permutation set 
for pcrmmomenta = l:3 
% clear data 
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NodeNumber = 1 ; 

clear internallines; i n t e r n al 1 i n e s ( NodeNumberMax ) =0; 
clear InternalTest ; In t er nalTes t =0; 
ltest =0; 

%% this evaluates all the graphs in the nodal complex 
for MontcRun = l:MonteRunSet ; 
for pcrmno = l:NodeNumberPerm 

%% j ranges from 1 to NodeNumber 

%% 

clear i n t e r n a 1 1 i n e s v al ; 

intcrnallinesval ( InternalLinesMax , d , NodeNumberMax , NodcNumbcrPcrm ) =0; 

Rainbowlntegral ( NodeNumberPerm ) =1; 

leftno =0; 

rightno =0; 

leftnodeno =0; 

rightnodeno =0; 

maxrightno =1; 

clear k ; 

%% store InternalLinesVal( nodenumber , linenumber , dimension ,permconfig) 
%% for all permutation sets contributing to the rainbow graphs 

for j =1: NodeNumberPerm 

clear momcntumsum ; 

momentumsum (NodeNumberMax , d) =0; 

clear NodeNumberMaxTwo=0; 
NodeNumberMaxTwo= ; 

for Dimension = l:d 

for NodcNumbcrTwo = 1 : NodeNumberMax 

for LineNumber = l: lineconfig (NodcNumbcrTwo , j ) 

if 1 i n e c o n f i g ( NodeNumbcrTwo , j ) >0 
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NodeNumberMaxTwo=NodeNumberTwo ; 
LineNumberTwo=LincNumber ; 

intcrnallinesval ( LineNumber , Dimension , NodcNumbcrTwo , j ) = (rand* 
kmax— kmax / 2 ) ; 

momontumsum ( NodcNumbcrTwo , Dimension ) =momcntumsum (NodeNumberT wo , Dimension) 
+ intcrnallinesval (LineNumber , Dimension , NodcNumbcrTwo , j ) ; 

% test = internallinesval( LineNumber , Dimension , NodeNumber , j 

); 

% propagatorsval( LineNumber , Dimension , NodeNumber ,j)=l/(tcst "2+m 

~2); 

end 

end 

end 
end 

end 

for NodeNumber = 1 : NodeNumbcrMax 

clear diff; diff(d)=0; 

if NodeNumber>l 

i f NodeNumber<NodeNumberMax+l 

for Dimension = l:d 

d i f f ( Dimension )=momentumsum ( NodeNumber — 1, Dimension) — 
momcntumsum ( NodeNumber , Dimension) ; 

momcntumsum ( NodeNumber , Di mens ion )=momcntumsum ( NodeNumber , 
Dimcnsion) + d i f f ( Dimension ) ; 

end 

for Dimension = l:d 

internallinesval ( LineNumberTwo , Dimension , NodeNumber , j )=d i f f ( 
Dimension) + intcrnallincsval ( LineNumberTwo , Dimension , 
NodeNumber , j ) ; 

end 
end 

end 

i f NodeNumber==l 
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for s = l:ppoint 

for Dimension = l:d 

internallinesval ( LineNumberTwo , Dimension , NodeNumber , j )= 

intcrnallincsval ( LineNumberTwo , Dimension , NodeNumber , j )+ 
kmomol ( s , Dimension , samplex ) ; 

momcntumsum ( NodeNumber , Dimension)=momentumsum ( NodeNumber , 
Dimension )+kmomol ( s , Dimension , samplex ) ; 

end 

end 

end 



% i f NodeNumber=NodcNumberMax 

% for s = l:qpoint 

% for Dimension = l:d 

% internallinesval( LineNumberTwo , Dimension , NodeNumber , j )= 

internallinesval ( NodeNumber , Dimension )+kmomor (s , Dimension , samplex) ; 

% momcntumsum ( NodeNumber , Dimension )=momcntumsum ( NodeNumber , 

Dimension )+kmomor (s, Dimension, samplex) ; 

% end 

% end 

% for Dimension=l:d 

% d i f f ( D i mens ion )=— momcntumsum ( NodeNumber — 1 , Dimension ) + 

momentumsum ( NodeNumber , Dimension)— internallinesval ( LineNumberTwo , 
Dimension , NodeNumber , j ) ; 

% momcntumsum ( NodeNumber , Dimension)=momentumsum ( NodeNumber , 

Dimension)+diff (Dimension) — internallinesval ( LineNumberTwo , Dimension , 
NodeNumber , j ) ; 

% end 

% for Dimension=l:d 

% internallinesval ( LineNumberTwo , Dimension , NodeNumber , j )=d iff ( 

Dimension) ; 

% end 

% end 
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%% end NodcNumbcrPerm 



end 

%% data of value of internal lines is stored in i n t er n al 1 i n e s v al 

%% the value of the propagators is stored in the propagators val 

%% momentum conservation is chosen to change the last entry which isnt 
implemented yet 

%% 

%% tree evaluation 

%% 

%% use a calling function 

%% 

%% at a node there are two sets of momenta 

%% internallinesval ( LineNumber , Dimension , NodeNumber , j ) 

%% internallinesval ( LineNumber , Dimension , NodeNumber + 1 , j ) 

%% care must be taken to define the external momenta of the rainbow gra 

%% as momenta is defined from the left to the node 



RainbowComputetwo=0; 



%% permutation set for the four point function 
i f permmomenta==l 

kmomol ; kmomor ; 

end 

%% l<->4 

i f permmomcnta==2 

testvec (1 ,: ,:) =kmomol ( 1 , : , : ) ; 
testvec (2 ,: ,:) =kmomor ( 2 , : , : ) ; 

kmomol (1 , : , : )=tcstvec ( 2 , : , : ) ; 
kmomor (2 ,: ,:)=tcstvcc (1 ,: ,:) ; 

end 

%% 4<->3 which is an overall l<->3 

i f permmomenta==3 

testvec (1 ,: ,:) =kmomol ( 1 , : , : ) ; 
testvec (2 ,: ,:) =kmomor ( 2 , : , : ) ; 

kmomol ( 2 ,: ,:)=tcstvcc (1 ,:,:); 



10 



kmomor (1 , : , : )=testvec ( 2 , : , : ) ; 

end 



% gnumber=8; 
% ppoint=2; 
% qpoint=2; 

% select the left and right momenta from a set kmomol and kmomor 

clear kmomoleft ; 
clear kmomoright ; 

kmomoleft ( ppoint ,d)=0; 
kmomoright ( qpoint , d) =0; 

%%%%% this is a subroutine to compute the rainbow iterations %%%%%% 
%% 

%% input external lines from the left side number p 
%% right side number q 

%% 

%% with the number of coupling constants 

%% this is the total number of internal lines 

%% 

%% the output is the value of the graphs subject to the number of internal 
lines 

%% 

%% summed over all their contributions 

%% 
% 

%%%%% these are the inputs 

%% partition into all possible rainbow graphs 
%% npoint tree has gnumber=npoint — 2 

%% rainbow graph with gnumber can have ppoint+2+qpoint+2+nodenumber*2 



%% max number of couplings 

%% with all 4— point couplings in the internal nodes 

%% with gnumber can have ppoint+qpoint + (gnumber— ppoint— qppoint 
-4) 

%% with only two nodes and no internal nodes 

%% 



% Two Nodes: InternalLines+ppoint — 2+InternalLines+qpoint — 2=gnumber 

InternalLinesMax = (gnumber+4— ppoint— qpoint ) /2; 

InternalLinesMin = 2; 
% NodeNumber is nodes minus one 

%% this is the max number of nodes all internal line numbers arc two 
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%% NodeNumbcrMax = (gnumber— ppoint— qpoint +2) /2 
NodeNumberMax = (gnumber— ppoint — qpoint +2) / 2 ; 
RainbowTree = 1.0; 
RainbowDcnom = 1.0; 
Atree = 1.0; 
Denominator = 1.0; 

for NodcNumbcr = 1 : NodcNumberMax — 1 

RainbowTree = 1.0; 
RainbowDcnom =1.0; 

leftno = lineconfig ( NodeNumber , permno ) ; 
rightno=lineconfig ( NodcNumbcr + 1 , permno ) ; 
if right n o >0 
if lcftno>0 

%% store in the number of lines to the right most node 
maxrightno=lineconfig ( NodeNumber + 1 , permno ) ; 

%% 

clear k; k ( 1 c f t n o+r ight no , d) =0; 

for LineNumber=l: leftno 
for dindcx = l:d 

k ( LineNumber , dindex)=internallinesval ( LincNumbcr , dindex , NodcNumbcr , 
permno ) ; 

end 

end 

for LincNumber = l: lincconfig ( NodcNumbcr + 1 , permno ) 
for dindex = l:d 

k ( LineNumber+lef t no , dindex ) = i n t e rn al 1 i n es v al (LineNumber , dindex , 
NodcNumbcr + 1, permno) ; 
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end 

end 

%%% this has setup the momenta variables in an array k_l k_2 etc k_n 
%%% for the tree calling function 

%%% k.l k.2 ... k.a left k.{a+l} ... k.n right 

%% call in the numerical values of the trees 

Atree=treecall(leftno+rightno ,k,N,truenumbertree) ; 
RainbowTree=RainbowTree * Atree ; 
Denominator=denomcall ( leftno+rightno ,k) ; 
RainbowDcnom=RainbowDenom * Denominator ; 
RainbowTree=RainbowTree *kmeasure " ( leftno —1) ; 

end 

end 

%% end the NodcNumbcr loop 
end 

%% the two external trees have to be included 
% first the left one 

lef t no=ppoint ; 

rightno=lincconfig (1 ,pcrmno) ; 
rightnodeno =1; 

clear k; k ( le f t n o+r ight no , d ) =0; 

for i =1: lef t no 

for dindex = l:d 

k ( i , dindcx )=kmomol ( i , dindcx , samplex ) ; 

end 

end 

for LineNumber=l: lineconfig (1 , permno ) 
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for dindex = l:d 

k ( LineNumber+lef t no , dindex) = internallinesval ( LineNumber ,dindcx, right node no , 
permno ) ; 

end 

end 

%% call in the numerical values of the left tree 

Atree=t r e e c a 1 1 ( 1 e f t no+right no ,k,N, truenumbertree) ; 

RainbowTree=RainbowTree * Atree ; 

Denominator=denomcall ( leftno+rightno ,k) ; 

RainbowDenom=RainbowDcnom * Denominator ; 
% the right tree 

rightno=qpoint ; 

for i = 1 : NodeNumberMax 

if lineconfig ( i , permno) >0; 

leftno=lineconfig(i , permno) ; 
leftnodeno = i ; 

end 

end 

clear k; k ( 1 c f t n o+r ight no , d) =0; 

for LincNumber=l: leftno 
for dindcx = l:d 

k ( LineNumber , dindcx)=internallinesval ( LineNumber ,dindex , leftnodeno , permno 

); 

end 

end 

for i=l:qpoint 

for dindex = l:d 

k(i + lcftno ,dindcx )=kmomor ( i , dindcx , samplex) ; 
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end 

end 

%% call in the numerical values of the right tree 

Atree=treecall ( leftno+rightno ,k,N, truenumbertree) ; 

RainbowTree=RainbowTree * Atree ; 

Denominator=denomcall ( leftno+rightno , k) ; 

RainbowDcnom=RainbowDcnom * Denominator ; 

RainbowTree=RainbowTree *kmeasure " ( leftno —1) ; 

%% this completes the tree nodal complex evaluation for the perm 
configuration 

%% have to sum over the individual graphs in the nodal complex this is 
%% denoted by permno 

% evaluate the graph 

% amputate the external lines 

dcnom = l; 

for i=l:ppoint 

dcnom = l/(kmomol ( i ,: , samplex ) *kmomol ( i ,: , samplex ) '+m" 2) *denom ; 

end 

for i=l:qpoint 

dcnom = l/(kmomor ( i ,: , samplex ) *kmomor ( i ,: , samplex ) '+m" 2) *denom ; 

end 

Rainbowlntegral ( permno )=RainbowTree *RainbowDcnom " (1/2) / denom " ( 1 / 2) ; 

% there is a z_2 sign ambiguity that needs to be determined from the 
% doubling of thepropagators in minkowski space 

%% end of permutation set 

end 

%% add up all the contributions from the multiple configurations 
Rainbow!ntegralScalar=0; 



15 



for j = 1 : NodeNumberPcrm 

Rai nbo win tegralScalar=Rainbo win tegralScalar+Rainbow Integral ( 

end 

RainbowComputetwo=RainbowIntegralScalar+RainbowComputetwo ; 

%% end of montc carlo simulation 
end 

%% end the permmomcnta run 
end 

RainbowCompute( samplex )=RainbowComputctwo ; 
%% end the samplex momenta run 
end 




%% t r e e c a 1 1 . m 

%% 
%% 

function [ treecall]=treecall (n , k ,N, truenumbertree) 
npoint=n ; 

%% 

%% numerical values of trees 

%%% 

%% external momenta are inputs k(j ,dindex) 

%% 

%% npoint is the number of legs j=l,...,n 

% 

%% N is the t_i~{[q]} matrix from which the trees are made 

%% 

»»K defaults 
m= 1 . ; 
d = 4; 

TotTreeValue = 0; 

clear kmominit ; kmominit (n , d) =0; 
clear kmomo; kmomo(n , d) =0; 
clear kindex ; kindex(d)=0; 
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pole =0; 

<mmmmmmmmmmmmmmmmmmi 

for a = l:npoint 
for j = 1 : d 

kmominit (a , j )=k ( a , j ) ; 
kmomo ( a , j ) =k ( a , j ) ; 
kindex ( j ) =0; 
end 
end 

% input a Nmatrix and output a number 
% sum over the trees 
for j =1 : tr uenumbertree ( npoint ) 
kprod = l; 
for a = l:npoint 

for b = l:npoint— 1 
for l=l:npoint 
for j =l:d 

kmominit (1 , j )=k(l , j ) ; 
kmomo( 1 , j )=k( 1 , j ) ; 
kindex ( j ) =0; 
end 
end 

pole=N( a ,b,j , npoint) ; 
if pole>0 

kindex ( : )=kmomo(a , : ) ; 
for c = l:b 

if a+c— Knpoint+1 
for e = l:d 
kindex ( c ) =kmomo ( a+c — 1 , e )+kindex ( c ) ; 
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end 

end 

if a+c— l>npoint 
for e = l:d 

kindex ( e )=kmomo( a+c— 1— npoint , e )+kindex ( c ) ; 
end 

end 

end 

kprod=kprod / (kindex* kindex '+m~ 2) ; 
end 

% storage of denominators as a function of steps 
end 
end 

TreeValue=kprod ; 
% add up the tree values 

TotTrccValuc=TotTreeValue+Tree Value ; 

end 

%% result for the treevalues is TotTreeValue 
%% rename to treecall 
treecall=TotTreeValue/2; 




%% NodalComplex .m 

%% 
%% 



function [ NodalComplex , NodcNumberPerm] = NodalComplex ( gnumber ,ppoint,qpoint) 

% clear data 
NodeNumber = 1 ; 
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NodeNumberMax = (gnumber— ppoint — qpoint +2) / 2 ; 
clear lineconfig ; clear lineconfigtwo ; 
li ncconf ig ( NodeNumberMax , 1 ) =0; 
lineconfigtwo ( NodeNumberMax , 1 ) = ; 
NodcNumbcrPcrm = 1 ; 
ltest =0; 
ptwo=0; 

kSampleMax = 400; 

% Two Nodes: I nt er nalL i nes+ppoint — 2+1 nt er nalL i nes +qpoint — 2=gnumbcr 
InternalLinesMax = (gnumber+4— ppoint— qpoint ) /2; 
InternalLinesMin = 2; 

%% configs are stored in the array internallines ( j ) 
for NodcNumber = l:NodeNumberMax 
for ktwo = l:kSampleMax 
clear internallines; 
for a = l:NodcNumbcr 

internallines ( a ) =0; 

end 

%% a random sampling to determine the configs 
for j =l:NodeNumber 

InternalTest=0; 
while InternalTest <2; 
InternalTest=round ( rand* InternalLinesMax) ; 

internallines (j )=InternalTcst ; 
end 

end 

addinternal =0; 
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for j = 1 : NodeNumber 

addinternal=addinternal+intcrnallincs ( j ) ; 

end 

%% relation between gnumber and node line configuration 

%% 

%% ppoint+internallinc 1 — 2 + i n t e r n al li n e 1 + i n t e r na 1 1 i ne 2 — 2 + ldots + 
%% intcrnalline (NodeNumber — l)+internallineNodeNumber —2 

%% + internallineNodeNumber+qpoint — 2 

%% 

%% gnumber = ppoint+qpoint —(NodeNumber + 1) *2+\sum" { NodeNumber } 2* 
internallinej 

%% 

%% all double internal gnumber=ppoint+qpoint — (NodeNumber + 1) *2+NodeNumber 
*2*2 

%% NodeNumber = (gnumber— ppoint— qpoint +2) /2 

addinternal =2* addinternal+ppoint+qpoint — 2*(NodeNumber + l) ; 
if addinternal=gnumber 
ptwo=ptwo + l; 

for a=l: NodeNumber 

lineconfigtwo ( a , ptwo ) = i n t c r n al 1 i n e s (a) ; 
end 

end 

%% end random sample 
end 

%% end NodeNumber loop 
end 

% 190 to 312 

%% eliminate the duplicates and store in lincconfig 
clear lincconfig; 
lincconfig ( NodeNumberMax , 1 ) = ; 
clear flagnumber ; clear flagnumbertwo ; 
flagnumber (ptwo) =0; flagnumbertwo (ptwo) =0; 
NodeNumberPerm = 1 ; 
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for r = l:ptwo — 1 
for s=r : ptwo 

tcstt wo =0; 

for j = 1 : NodcNumberMax 

testt wo=testtwo+abs (lineconfigtwo (j , r ) — lineconfigtwo (j , s 

end 

if testtwo==0 

% if r<s 

if flagnumbertwo ( s )==0 

flagnumbcr(r)=l; 
flagnumbertwo ( s ) =1; 

end 

% end 
end 

end 

end 

for r = l:ptwo 

if flagnumber ( r )==1 

lincconfig (: , NodcNumberPerm ) = lineconfigtwo (: ,r) 
NodeNumberPerm=NodcNumbcrPcrm + 1 ; 

end 

end 

NodcNumbcrPcrm=NodeNumberPerm — 1 ; 

%%% this completes the possible permutations of nodes and internal lines 
%%% this order and NodcNumberPerm is the number of permutations 

%%% all possible node configurations and internal lines arc stored in 
%% the variable 

%%% 

%%% lineconfig (NodeNumber , a) 

%%% 

%%% NodeNumber=node number 

%%% 

%%% a=configuration number up to NodcNumberPerm 

%% 
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%%% 



NodalComplex=l i n e c o nf i g ; 
NodalComplcxPcrm=NodcNumbcrPcrm ; 



%% 
%% 

%% NMatrix.m 

%% 
%% 



function [ NMatrix , t ruenumbert ree] = NMatrix (npoint , npoint lower) 
n=npoint ; 

% npointlower is another useful delimiter 

%% 

%% t_i"{[q]} formation for the matrix of trees 

%% 

%% npoint is the required dimension and is input 

%% 

%% matrix N is spanned by N(i ,q,j ,m) with m corresponding to m-point 
trees 

%% 

%% i,q represent t _ i " { [ q ] } 

%% 

%% j corresponds to the different trees up to truenumbertree (m) 

%% 



clear N; 

clear treenumber ; 
treenumber ( npoint ) =0; 

for spoint=npointlower : npoint 

rpoint=npoint — spoint+npoint lower ; 
n=r point ; 



clear sigmanumber ; clear flagnumber ; 



for k = l :Gamma( n+1) 

sigmanumber ( 1 , k )=n ; 
for j=2:n-2 
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sigmanumbcr ( j ,k)=0; 

x = l; 
while x<2 

x=round (n*rand) ; 

sigmanumbcr ( j ,k)=x; 
end 

end 
end 

% reorder the entries and eliminate redundancy 

for k = l:Gamma(n+l) 
for m=l :Gamma( n+1) 
for j = l:n— 3 

if sigmanumbcr ( j , k)<sigmanumbcr ( j +1 ,k ) 

test l=sigmanumbcr ( j , k) ; t est 2=sigmanumber ( j +1 ,k) ; 
sigmanumber (j , k)=test 2 ; sigmanumbcr (j+l,k)=testl ; 
end 
end 
end 
end 

%% the redundancy in sigmanumber is eliminated 
%% values are stored in rp 
%% note that the zero vector is allowed 
v=size (sigmanumbcr ) ; dimp=v(2); 

for i=l:dimp 
flagnumber ( i ) =0; 

end 
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for i=l:dimp 

% redundancy check 
for m=l:dimp 

i f m=i 

flagnumber (m) =1; 

end 



i f m~= i 

i f sigmanumber ( : ,m)=sigmanumbcr ( : , i ) 
if i >m 

flagnumber (m) =0; 

else 

flagnumber (m) =1; 

end 
end 

end 

end 
end 

totalflagno=0; 
for i=l:dimp 

testmin =0; 
for j=l:n-2 
if sigmanumber (j , i )==2 

testmin=testmin + 1; 
end 
end 

if testmin>l 
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flagnumber( i)=0; 

end 

totalflagno=totalflagno+flagnumber ( i ) ; 
end 

rp ( 1 , 1 ) =0; clear rp 

for i =1: totalflagno 

for j =l:n— 2 
rp(j , i)=0; 

end 
end 

treenumbertwo =0; 
k=0; 
for i=l:dimp 

% eliminate duplicate graphs 

if flagnumbcr ( i )==1 
k=k+l; 

rp ( : , k)=sigmanumber ( : , i ) ; 
treenumbertwo=treenumbertwo + 1; 
end 
end 
P=rp; 

treenumber ( rpoint )=treenumbertwo ; 

nmmmmmmmmm^mmwwwo 
%% 

Wo 2) deduce the t_i"{[q]} variables 

Wo 

%% store in matrix T( i , q , treeno ) 

%% 

clear P; 
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clear M; 

M(npoint , npoint , treenumber ( npoint ) ) =0; 

for trecno =1 : treenumber ( rpoint ) 
for i =l:n 

for q=l:n— 1 

M( i , q , treeno * rpoint ) =0; 
end 
end 
end 



for treeno =1 : treenumber ( rpoint ) 
% for b=0:n-2 
b = 0; 

%% first construct (p , [ p_m] ) sets, in P(i,j, treeno) matrix 

for j =l:n 

P(j ,trceno)=0; 
Qo (j ,treeno)=0; 
PP(j , trccno)=0; 
Q( j , treeno ) =0; 

end 

for j =l:n 

for i=l:n-2 

if rp ( i , trecno )=j 

P(j ,treeno)=P(j ,trecno)+l; 
Qo ( j , t r e e n o ) =Qo (j ,treeno)+l; 
Q(j , treeno )=Qo( j , treeno ) ; 
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end 

if j<n 

if P(j , treeno)>0 
for k= j + 1 : n 

PP ( k , t r e e n o ) =P (j , tree no); 
end 
end 
end 

end 
end 

for j =l:n 

if P(j, treeno) ==0 

P ( j , treeno) =PP ( j , t r e e n o ) ; 

end 
end 



%% the numbers labeling the graphs are stored in rp(i,j) with first entry n 
%% T is initialized to zero 

%% to find a pole n— l=number of nodes : search for subtrees 
%% scan all initial lines 

%% f i n a 1 1 i n e — i n i t i a 1 i n e — l=nodes in tree between finalline and initialline 
m=n ; 

nodesl=0; nodes2=0; 
b = 0; 

% permute the roundabout labeling as 1— >permn 

for j = 1 : n 

i f j +b<n+l 
Q(j , treeno )=Qo( j+b , treeno ) ; 
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end 

i f j +b>n 
Q(j , treeno )=Qo( j+b— n , trccno ) ; 
end 

end 

for i n i t i a 1 1 i n e =1 : n— 1 

for finalline = initialline+l:n 

nodesl=0; nodes2=0; 

init count=Q( i nit i alii ne , treeno); 
finalcount=Q( final line , treeno); 

Qsum=0; 

% around clockwise the tree diagram 
if f in alii ne — i n i t i a 1 1 i n e >1 
if initcount ==0 
if finalcount >0 

for k = l: finalcount 
Qsum=0; 

for alpha= initialline+l:finalline — 1 
Qsum=Qsum+Q ( alpha , treeno ) ; 

end 

Qsum=Qsum+k ; 

% skipping rule applied here 

if Q(finalline — 1, treeno )==0 
% check for subtree 

if Qsum=f i n a 1 1 i n c — i n i t i a 1 1 i n e 
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i=i n i t i a 1 1 i n e ; 

q=f i nalline— initialline+1; 

M(i,q,treeno) =b + 1 ; 

end 

end 
end 

end 
end 
end 

% f i n alii n e — i n i t i a 1 1 i n e =1 case 

if finalline— initialline ==1 
if Q( initialline , trccno) ==0 
if Q(finalline ,treeno)>0 
i= initialline ; 
q= finalline — initialline +1; 
M(i ,q,treeno) =b + 1 ; 

end 

end 
end 

% end of initalline and finalline loop 
end 
end 

% end of the permutation 
% end 

% end the tree graph count 
end 

% no end to the spoint loop yet 
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% permute the matrices 
% wasteful of memory 
for i =1: rpoint 

treeno=treenumber ( rpoint ) ; 

for s = l:trcenumber ( rpoint ) 
for a = l:rpoint 



for b = l:rpoint 

if a+i — Krpoint +1 

M(a+i -l,b, ( i - l)*trccno+s )=M(a,b, s ) ; 

end 

if a+i— 1> rpoint 

M(a+i —1— rpoint ,b,(i— l)*trceno+s )=M( a , b , s ) ; 

end 



end 
end 
end 

end 

% for tr ccno =1 : rpoint * treenumber ( rpoint ) 
% for i=l:n 
% for q=l:n-2 
% if M(i ,q,treeno )=1 

% if i+q<n 

% M( i+q+l,n—q , treeno ) =0; 



% 



end 



% 



if i+q>n 



% 



M( i+q— n , n— q , treeno ) =0; 



% 



end 



% 



end 
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% end 

% M(i ,rpoint ,treeno)=0; M( i , 1 , treeno ) =0; M( i ,rpoint — l,trceno)=0; 
% end 
% end 

for j=l:rpoint*treenumbcr (rpoint) 

for k=j + l:rpoint*treenumber (rpoint ) 

test =0; 
for s = l:rpoint 

for t = l:rpoint 

test=abs (M(s , t , j )-M(s , t , k) )+tcst ; 

end 

end 

if test==0 
if k>j 

M(l ,1 ,k)=M(l ,1 ,k)+43*rand; 
end 

end 

end 
end 

% there could be an error in the t_i"{[q]} matrix but not likely 

Nbefore ( rpoint , rpoint , rpoint *treenumber ( rpoint )) =0; 
N(rpoint , rpoint ,rpoint*treenumbcr(rpoint) ,rpoint)=0; 
nindex = 0; 

for s = l:rpoint*treenumber (rpoint) 
sumM=0; 
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for i =1: rpoint 
for j =1: rpoint 

sumM=sumM4M( i , j , s ) ; 

end 
end 

if sumM=rpoint — 3 

nindex=nindex + l; 
Nbeforc(: ,: ,nindex )=M( : , : , s ) ; 

end 

N(: ,: ,: , rpoint )=N before (:,:,:); 
truenumbertree( rpoint )=nindcx ; 
end 

%%% end loop on the rpoint index 
end 

mmmmmmmmmmmmmmmmmmmmm 
%% 

%% innerproduct .m 

%% 
%% 

function [innerproduct] = innerproduct (ppoint ,qpoint , kmomoleft , kmomoright ) 
load Nmatrix 

for j=l:ppoint 

test=test+kmomoleft ( j , : ) *kmomoright ( j ,:) '; 

end 

innerproduct=test ; 

mmmmmmmmmmmmmmmmmmmmi 
%% 

%% denomcall .m 

%% 
%% 
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function [denomcall] = denomcall (numberlines , k) 



m=1.0; 
dcnom = 1 ; 
d = 4; 

for i=l:d 
p(i)=0; 
end 

for j =1: numberlines 
for i=l:d 

q(i )= k (j , i ) ; 

end 

denom=denom / ( q*q'+m" 2) ; 

end 

denomcall=denom ; 
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